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Abstract 

We  present  a  simple  0(nm  +  n-^  log  U)  sequential  algorithm  for  the 
maximum  flow  problem  with  n  nodes,  m  arcs,  and  a  capacity  bound  of  U 
among  arcs  directed  from  the  source  node.  Under  the  practical  assumption 
that  U  is  polynomially  bounded  in  n  ,  our  algorithm  runs  in  time  0(nm  + 
n-^  log  n).  This  result  improves  the  current  best  bound  of  0(nm  log  (n^/m), 
obtained  by  Goldberg  and  Tarjan  ,  by  a  factor  of  log  n  for  non-sparse  and  non- 
dei  se  networks  without  using  any  complex  data  structures. 
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The  maximum  flow  problem  is  one  of  the  most  fundamental  problems 
in  network  flow  theory  and  has  been  investigated  extensively.   This  problem 
was  first  formulated  by  Ford  and  Fulkerson  [1956]   who  also  solved  it  using 
their  well-known  augmenting  path  algorithm.    Since  then  a  number  of 
algorithms  have  been  developed  for  this  problem,   as  tabulated  below.   Here  n 
is  the  number  of  nodes,  m  is  the  number  of  arcs,  and  U  is  an  upper  bound 
on  the  integral  arc  capacities. 

#  Due  to  Running  Time 

1  Ford  and  Fulkerson  [  1 956]  0(nm  U) 

2  Edmonds  and  Karp  [1972]  0(nm2) 

3  Dinic  [1970]  O(n^) 

4  Karzanov  [1974]  OCn^) 

5  Cherkaskyll977]  ©(n^m^/^) 

6  Malhotra,  Kumar  and  Maheshwari  [1978]  O(n^) 

7  Galil[1978]  0(n5/3m2.'3) 

8  Galil  and  Naamad  [1978]  0(nm  log^  n) 

9  Shiloach  and  Vishkin    [1982]  OCn^) 

10  Sleator  and  Taijan  [1983]  0(nm  log  n) 

11  Tarjan[1984]  OCn-"*) 

12  Gabow[1985]  O(nmlogU) 

13  Goldberg  [1985]  OCn^) 

14  Goldberg  and  Tarjan  [1986]  0(nm  log  (n^/m)) 


Edmonds  and  Karp  [1972]  showed  that  the  Ford  and  Fulkerson  [1956] 
algorithm  runs  in  time  0(nm-^)  if  flows  are  augmented  along  shortest  paths 
from  source  to  sink.   Independently,  Dinic  [1970]  introduced  the  concept  of 
shortest  path  networks,  called   layered  networks,  and  obtained  an  O(n-m) 
algorithm.   This  bound  was  improved  to  O(n^)  by  Karzanov  [1974]   who 


introduced  the  concept  of  preflows   in  a  layered  network.    A  preflow  is  similar 
to  a  flow  except  that  the  amount  flowing  into  a  node  may  exceed  the  amount 
flowing  out  of  a  node.    Since  then,  researchers  have  improved  the  complexity 
of  Dinic's  algorithm  for  sparse  networks  by  devising  sophisticated  data 
structures.    Among  these  contributions,   Sleator  and  Tarjan's  [1983]  dynamic 
tree  is  the  most  attractive  from  a  worst  case  point  of  view. 

The  algorithms  of  Goldberg  [1985]  and  of  Goldberg  and  Tarjan  [1986] 
were  a  novel  departure  from  these  approaches  in  the  sense  that  they  did  not 
construct  layered  networks.     It  contains  the  essence  of  Karzanov's  preflow 
method  but  does  not  maintain  layered  networks.  Their  algorithm  maintains 
a  preflow  and  proceeds  by  pushing  flows  to  nodes  estimated  to  be  closer  to  the 
sink.     To  estimate  which  nodes  are  closer  to  the  sink,  it  maintains  a  distance 
label  for  each  node  that  is  a  lower  bound  on  the  length  of  a  shortest 
augmenting  path  to  the  sink.  Distance  labels  are  a  better  computational  device 
than  layered  networks  since  these  labels  are  simpler  to  understand,  easier  to 
manipulate,   and  easier  to  use  in  a  parallel  algorithm.     Moreover,  by  cleverly 
implementing  the  dynamic  tree  data  structure,  Goldberg  and  Tarjan  obtained 
the  best  computational  complexity  for  sparse  as  well  as  dense  networks. 
Bertsekas  [1986]    recently  developed  an  algorithm  for  the  minimum  cost  flow 
problem  that  simultaneously  generalizes  both  the  Goldberg-Tarjan  algorithm 
and  Bertsekas  [1979]  auction  algorithm  for  the  assignment  problem. 

For  problems  with  arc  capacities  polynomially  bounded  in  n  ,  our 
maximum  flow  algorithm  is  an  improvement  of  Goldberg  and  Tarjan's 
algorithm  and  uses  concepts  of  scaling  introduced  bv  Edmonds  and  Karp 
[1972]  for  the  minimum  cost  flow  problem  and  later  extended  by  Gabow  [1985] 


for  other  network  optimization  problems.    The  bottleneck  operation  in  the 
straightforward  implementation  of  Goldberg  and  Tarjan's  algorithm  is  the 
number  of  "non-saturating  pushes"  which  is   0(n-^)  .    However,  they  reduce 
the  computational  time  to  0(nm  log  (n^/m))  by  a  very  clever  application  of 
the  dynamic  tree  data  structure.   We  show   that  the  number  of  non-saturating 
pushes  can  be  reduced  to    0(n^  log  U)  by  using  "excess  scaling."  Our 
algorithm  performs   log  U   scaling  iterations;   each  scaling  iteration  requires 
O(n^)  non-saturating  pushes  if  we  push  flows  from  nodes  with  sufficiently 
large  excesses  to  nodes  with  sufficiently  small  excesses    while  never  allowing 
the  excesses  to  become  too  large.      Consequently,   the  computational  time  of 
our  algorithm  is   0(nm  +  n^  log  U).   Under  the  reasonable  assumption  that 
U  =  OCn*^^^^)     (i.e.,   it  is  polynomial  in  n)  ,  the  algorithm  runs  in  time 
0(nm  +  rr  log  n)   which  improves  the  bound  of  Goldberg  and  Tarjan's 
algorithm  by  a  factor  of  log  n   for  both  non-sparse  and  non-dense  networks, 
i.e.,  networks  for  which  m  =  OCn^-"^^)    and    m=  Q  (n^"^^)   for  some    e     with 
0  <  E  <  1    .     Moreover,  our  algorithm  is  easy  to  implement  (and  is  much 
more  efficient  in  practice),  since  it  requires  only  elementary  data  structures 
with  little  computational  overhead. 

Our  algorithm  is  computationally  attractive  even  if  U  is  not   0(n'^'^^) 
and  the  arc  capacities  are  exponentially  large  numbers.   In  this  case,  the 
uniform  model  of  computation,  in  which  all  arithmetic  operations  take   0(1) 
steps,  is  arguablv  inappropriate.   It  is  more  realistic  to  adopt  the  logarithmic 
model  of  computation  (as  described  by  Cook  and  Reckhow  [1973])  which 
counts  the  number  of   bit  operations.    In  this  model,  most  arithmetic 
operations  take  Odog  U)  steps  than  0(1)  steps. 

Using  the  logarithmic  model  of  computation  and  modifying  our 
algorithm  slightly  to  speed  up  arithmetic  operations  on  large  integers  ,   we 


obtain   an  0(nm  log  n   +  n-  log  n   log  U)  algorithm  for  the  maximum  flow- 
problem.    The  corresponding  time  bound  for  the  Goldberg-Tarjan  algorithm 
is    0(nm  log  (n^/m)  log  U).     Hence  as  U  becomes  exponentially  large,  our 
algorithm  surpa         T^oldberg  and  Tarjan's  algorithm  by  a  factor  of  m/n  in 
the  non-dense  c:  n  other  words,  the  relative  worst  case  performance  of 

our  algorithm  i?  "kingly  superior  as   U  increases  in  size.   Our  results  on 

this  logarithmic  moaei  of  computation  will  be  presented  in  Ahuja  and  Orlin 
[19S7a] . 

Our  maximum  flow  algorithm  is  difficult  to  make  "massively  parallel" 
since  we  push  flow  from  one  node  at  a  time.   Nevertheless,  with  d  ==  f  m/n1 
parallel  processors  we  can  obtain  an  0(n^  log  Ud)  time  bound.  Under  the 
assumption  that  U  =   ©(n"^'^^)  ,  the  algorithm  runs  in  0(n-   log  n)  time,  which 
is  comparable  to  the  best  available  time  bounds  obtained  by  Shiloach  and 
Vishkin  [19S2]    and  Goldberg    and  Tarjan  [1986]  using  ?:  ;7flrfl//f/  processors. 
Thus,    our  algorithm  has  an  advantage  in  situations  for  which  parallel 
processors  are  at  a  premium.   Our  work  on  the  parallel  algorithms  will  be 
presented  in  Ahuja  and  Orlin  [19S7b]  . 


1.        Notation 

Let  G  =  (N,  A)  be  a  directed  network  with  a  positive  integer  capacity  Ujj 

for  even,'  arc  (i,  j)  e  A.  Let  n  =  INI    and  m  =  I A  I  .  The  source  s     and    sink 
t      are  two  distinguished  nodes  of  the  network.   We  assume  without  loss  of 
generality  that  the  network  does  not  contain  multiple  arcs.    We  further 
assume  that  for  everv  arc  (i,  j)  e  A  ,  an  arc  (j,  i)  is  also  contained  in  A  , 

possibly  with  zero  capacity.  Let  U  =  max      (UgJ. 

(s,  j)6  A 


A  flow  is  a  function    x:  A  — >  R   satisfying 


S         x^i        -  X        Xjj    =  0  ,  forall    i  e  N  -  {s,  t},  (1) 

{j:(i,i)6A}  {j:(i,j)GA} 


I      x^t    =   V  ,  (2) 

{j:  (j,  t)  e  A) 


0  <  Xjj  <  u-j    ,       for  all     (i,  j)  e  A  ,  (3) 


for  some    v  >  0  .      The  maximum  flow  problem  is  to  determine  a  flow  x      for 
which    V     is  maximized. 

A  -prcfloiv     X    is  a  function    x:  A  — >  R    which  satisfies   (2)  ,  (3),   and 
the  following  relaxation  of   (1): 


I         Xji        -  I        Xjj    >  0    ,  for  all    i  e  N  -  {s,  t}.  (4) 

{j:  (j,i)€  A}  {j;  (i,j)e  A} 


The  algorithms  described  in  this  paper  maintain  a  preflow  at  each 
intermediate  stage. 

For  a  given  preflow  x  ,  we  define  for  each  node     i  g  N  -  (s,  t]  ,    the 
excess 


e,  =  I         Xjj        -  I       Xjj 

{j:  (j.i)e  A}  {j:  (i,j)e  A} 


A   node  with  positive  excess  is  referred  to  as  an   active  node  .     The 
residual  capacity  of  any  arc  (i,  j)  e  A    with  respect  to  a  given  preflow  x  is 
given  by      rj:  =         -  x^:    +  x-   .     The  network  consisting  only  of  arcs  with 

positive  residue  > cities  is  referred  to  as  the  residual  network.     Figure  1 

illustrates  thest  tions. 


Network  vsnth  arc  capacities. 

Node  1    is  the  source.    Node  4  is  the 

sink.  (Arcs  v»-ith  zero  capacities  are  not 

shewn) 


Network  with  a  preflow  x 


The  residual  network  with 
residual  arc  capacities 


Figure  1.       Illustrations  of  preflow  and  residual  network. 


A  distance  function    d:  N— >  Z"^   for  a  preflow  x    is  a 
function  from  the  set  of  nodes  to  the  non-negative  integers  .       We  say  tliat  the 
distance  function  is  valid  if  it  also  satisfies  the  following  two  conditions: 

CI.       d(t 

C2.       d(i;  ')  +  1  ,  for  every  arc  (i,  j)  e  A    with    r-  >  0  . 

Our  algorithm  maintains  a  vaUd  distance  function  at  each  iteration, 
and  labels  each  node  i   with  the  corresponding  value  d(i).    The  distance    d(i) 
is  a  lower-bound  on  the  length  of  any  path  from  node  i   to   t  in  the  residual 
network.     This  fact  is  easy  to  demonstrate  using  induction  in  the  value  d(i). 
If  for  each  i,  the  distance  label   d(i)   equals  the  minimum  length  of  any  path 
from   i  to  t    in  the    residual  netw^ork,   then  we  call  the  distance  label  exact. 
For  example,  in  Figure  1(c),  d  =  (0  ,  0  ,  0  ,  0  )    is  a  valid  distance  label,  though 
d  =  (3, 1,  2,  0)      represents  the  exact  distance  labels. 

We  define  the  arc  adjacency  list  A(i)  of  a  node    i  e  N    as  the  set  of  arcs 
directed  out  of  the  node  i  ,  i.e.,  A(i)  :  =   {(i,  k)  e  A  :  k  e  N}.    Note  that  our 
adjacency  list  is  a  set  of  arcs  rather  than  the  more  conventional  definition  of  a 
list  as  a  set  of  nodes. 

All  logarithms  in  this  paper  are  assumed  to  be  of  base  2  unless  stated 
otherwise. 
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2.        The  Preflow-Push  Algorithms 

The  preflozv-push  algorithms  for  the  maximum  flow  problem 
maintain  a  preflow  at  every  step  and  proceed  by  pushing  the  excess  of  nodes 
closer  to  the  sink.   The  first  preflow-push  algorithm  is  due  to  Karzanov  [1974]. 
Tarjan  [1984]   has  suggested  a  simplified  version  of  this  algorithm.   The  recent 
algorithms  of  Goldberg  [1985]  and  Goldberg  and  Tarjan  [1986]  are  based  on 
ideas  similar  to  those  presented  in  Tarjan  [1984]  but  use  distance  labels  to 
direct  flows  closer  to  the  sink  instead  of  constructing  layered  network.    We 
thus  call  their  algorithm  the  distance-directed  preflow-push  algorithm.    In 
this  section,  we  review  the  basic  features  of  this  algorithm,  which  for  the  sake 
of  brevity,  we  shall  simply  refer  to  as  the  preflow-push  algorithm.   Here  we 
describe  the   1-phase  version  of  the  preflow-push  algorithm  presented  by 
Goldberg  [1987].   The  results  in  this  section  are  due  to  Goldberg  and  Tarjan 
[1986]. 

All  operations  of  the  preflow-push  algorithm  are  performed  using 
only   local  information.    At  each  iteration  of  the  algorithm  (except  at  the 
initialization  and  at  the  termination)    the  network  contains  at  least  one  active 
node,  i.e.,  a  (non-source)  node  with  positive  excess.    The  goal  of  each 
iterative  step  is  to  choose  some  active  node  and  to  send  its  excess  "closer"  to 
the  sink,  with  closer  being  judged   with  respect  to  the  current  distance  labels. 
If  excess  flow  at  this  node  cannot  be  sent  to  nodes  with  smaller  distance  labels, 
then  the  method  increases  the  distance  label  of  the  node.    The  algorithm 
terminates  when  the  network  contains  no  active  nodes.    The  preflow-push 
algorithm  uses  the  following  subroutines: 
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PRE-PROCESS.  On  each  arc  (s,  j)  e  A(s)  ,  send  u^:  units  of  flow.  Let 
d(s)  =  n  and  d(t)  =  0  .  Let  d(i)  =  1  for  each  iit  s  or  t  . 
(Alternatively,  the  distance  label  for  each  node  i  ^^  s,  t  can  be 
detern^int  a  breadth  first  search  on  the  residual  network  starting  at 

node  t.  ) 

SELECT.         Select  an  active  node    i  ^^  s,  t    (i.e.,    Cj  >  0)  . 

PUSH(i).  Select  an  arc  (i,  j)  e  A(i)  with  r^j  >  0  and  d(i)  =  d(j)  +  1. 
Send   6  =   min  {e, ,  r^:  }  units  of  flow  from  node    i    to    j. 

We  say  that  a  push  of  flow  on  arc  (i,  j)    is  saturating  if    5  =  rjj  , 
and    non-saturating    otherwise. 

RELABEL(i).  Replace  d(i)  by  min{  dQ)  +1:  (i ,  j)  g  A(i)  and  r,:  >  0  }. 
Figure  2  contains  the  generic  version  of  the  preflow-push  algorithm. 


begin 

PREPROCESS 

while    there  is  an  active  node  in  N  -  {s,  t}    do 

begin 

SELECT   {let  i   denote  the  node  selected); 

end; 

If    there  is  an  arc  (i,  j)  €   A(i)    with    r^j  >  0  and 

d(i)  =  d(j)  +  l   then     PUSH(i) 
else  RELABEL(i) 
end; 

Figure  2.    The  preflow-push  algorithm 
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Figure  3  illustrates  the  steps    PUSH(i)     and     RELABEL(i)  as  applied  to 
the  network  in  Figure  1(a).   Figure  3(a)  specifies  the  preflow  determined  by 
PRE-PROCESS.    The  SELECT  step  selects  node  2  for  examination.  Since  arc 
(2,  4)  has  residual  capacity    r24     =    1    and        d(2)  =  d(4)  +  1  ,  the  algorithm 

performs   a  (saturating)  push   of  value   5  =  min{  1,  1}   units.  The  push 
reduces  the  excess  of  node  2  to  1  .  Arc  (2,  4)  is  deleted  from  the  residual 
network  and  arc  (4,  2)  is  added  to  the  residual  network.   Since  node  2  is  still 
an  active  node,  it  can  be  selected  again  for  further  pushes.  The  arcs  (2,  3)  and 
(2,  1)   have  positive  residual  capacities,  but  they  do  not  satisfy  the  distance 
condition.   Hence  the  algorithm  performs    RELABEL(2)  ,and  gives  node  2  a 
new  distance  d'(2)  =  min   {d(3)  +  1  ,  d(l)  +  1)  =  min  {2,  5}  =  2  . 

The  pre-process  step  accomplishes  several  important  tasks.   First,  it 
causes  the  nodes  adjacent  to   s   to  have  positive  excess,  so  that  we  can  select 
some  node  with  positive  excess.     (Otherwise,  the  algorithm  could  not  get 
started.)    Second,  by  saturating  arcs  incident  to  s  ,  the  feasibilit}'  of  setting 
d(s)  =  n   is  immediate.   Third,  since  d(s)  =  n      is  a  lower  bound  on  the  length 
of  the  minimum  path  from    s   to   t  ,   there  is  no  path  from   s   to   t  .      Since 
distances  in    d    are  non-decreasing  (see  Lemma  1   to  follow),  we  are  also 
guaranteed  that  in  subsequent  iterations  the  residual  network  will  never 
contain  a  directed  path  from   s  to  t,  and  so  there  can  never  be  any  need  to 
push  flow  from   s  again. 
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d(3)  =  l 


d(l)  =  4 


d(4)  =  0 


(a)      Thv  .'ork  with  a  preflow  after  the  pre-processing  step. 


d(3)  =  l 


d(l)  =  4 


d(4)  =  0 


(b)     After  the  execution  of  step  PUSH(2 


d(3)  =  1 


d(l)  =  4 


d(4)  =  0 


dt2)  =  2 


e      =  1 

-I 


(c)    After  the  execution  of  step  RELABEL(2). 
Fisurt  3.      Illustriitioiis  of  Push    and    Relabel    steps. 
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In  our  improvement  of  the  preflow-push  algorithm,   we  need  a  few  of 
the  results  given  in  Goldberg  and  Tarjan  [1986].   We  include  some  of  their 
proofs  in  order  to  make  this  presentation  more  self-contained. 

Lemma  1.      The  generic  preflow-push  algorithm  maintains  valid  distance 
labels  at  each  step.  Moreover,  at  each  relabel  step  the  distance  label  of  a  node 
strictly  increases. 

Proof.  First  note  that  the  pre-process  step  constructs  valid  distance 

labels.     Assume  inducti\  el}^  that  the  distance  function  is  valid  prior  to  an 
operation,  i.e.,  satisfies  the  validity  conditions  CI     and  C2  .   A  push 
operation  on  the  arc  (i,  j)  may  create  an  additional  arc   (j,  i)  with  r-  >  0  ,  and 

an  additional  condition       d(j)  <  d(i)  +  1      needs  to  be  satisfied.  This  validity 
conditions  remain  satisfied  since  d(i)  =  d(j)  +  1   by  the  property  of  the  push 
operation.   A  push  operation  on  arc   (i,  j)  might  delete  this  arc  from  the 
residual  network,  but  this   does  not  affect  the  validity  of  the  distance  function. 
During  a  relabel  step,  the  new  distance  label  of  node  i  is    d'(i)  =  min{d(j)  +  1: 
(i,  j)  e  A(i)  and    r-  >  0)  ,   which  is   consistent  with  the  validity  conditions. 

The  relabel  step  is  performed  when  there  is  no  arc  (i,  j)  €  A(i)    with 

d(i)  =  d(j)  +  1   and  r-:  >  0  .    Hence, 

d(i)  <   min{d(i)  +  l:    (i,  j)  e  A(i)  and    r- >  0}  =  d'(i)     ,  thereby  proving  the 

second  part  of  the  lemma.    ■ 


Lemma  2.  At  any  stage  of  the  preflow-push  algorithm,  each  node  with 
positive  excess  is  connected  to  node  s  by  a  directed  path  in  the  residual 
network. 
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Proof.  By  the  flow  decomposition  theory  of  Ford  and  Fulkerson  [1962]  , 

any  preflow  x  can  be  decomposed  with  respect  to  the  original  net^vork  G 
into  non-negative  flows  along    (i)  paths  from  source  s   to  t,   (ii)  paths  from  s 
to  active  nodes  pi'       and  (iii)  the  flows  around  directed  cycles.    Let  i  be  an 
active  node  in  th  How  x  of  G   .  Thus  there  must  be  a  path  P    from  s  to  i 

in  the  flow  decc  tion  of  x  .  Then  the  reversal  of  P  (P  with  the 

orientation  of  each  arc  reversed)   is  in  the  residual  network   G',   and   hence 
there  is  a  path  from    i   to  s  in  G'.    ■ 

Corollary  1.    For  each  node  i  e  N  ,    d(i)  <  2n. 

Proof.        The  last  time  that  node  i  was  relabeled,  it  had  a  positive  excess,  and 

hence  the  residual  network  contained  a  path  of  length  at  most   n-  1     from  i   to   s. 
The  fact  that  d^^  =  n  and  condition  C2  imply  that     dj  <  dg  +  n  -  1  <  2n.     ■ 

Corollary  2.      The  number  of  relabel  steps  is  less  than   2n'^  . 

Proof.  Each  relabel  step  increases  the  distance  label  of  a  node  by  at  least 

one,  and  bv  Corollarv  1  no  node  can  be  relabeled  more  than  at  most  2n  times.   ■ 

Corollarv  3.      The  number  of  saturating  pushes  is  no  more  than  nm  . 
Proof.  Suppose  that  arc  (i,  j)   becomes  saturated  at  some  iteration  (at 

which   d(i)  =  d(j)  +  1).   Then  no  more  flow  can  be  sent  on  (i,  j)   until  flow  is 
sent  back  from  j  to  i  (at  which  d'(j)  =  d'(i)  +  1   >  d(i)  +  1  =  d(j)  +  2);    This  flow 
change  cannot  occur  until   d(j)   increases  by  at  least  2.  Thus  by  Corollary  1 
arc  (i,  j)  can  become  saturated  at  most  n  times,  and  the  total  number  of  arc 
saturations  is  no  more  than   nm.     (Recall  that  we  assume  that    (i,  j)   and   (j,  i) 
are  both  in   A  ,   so  the  number  of  arcs  in  the  residual  network  is    no  more 
than  m  .  )      K 
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Lemma  3.      The  number  of   non-saturating  pushes  is  at  most     Zn-^m. 
Proof.  See  Goldberg  and  Tarjan  [1986]  . 

Lemma  4.      The  algorithm  terminates  with  a  maximum  flow. 
Proof.  When  the  algorithm  terminates,  each  node  in   N  -  {s,  t)  has  zero 

excess;    so  the  final  preflow  is  a  feasible  flow.  Further,  since  the  distance 
labels  satisfy  the  conditions  CI   and  C2    and    d(s)  =  n  ,  it  follows  upon 
termination  the  residual  network  contains  no  directed  path  from  s   to     t  . 
This  condition  is  the  classical  termination  criteria  for  the  maximum  flow 
algorithm  of  Ford  and  Fulkerson  [1962]  .      ■ 

The  potential  bottleneck  operation  in  the  preflow  based  algorithms  due 
to  Karzanov  [1974]  ,  Tarjan  [1984],  and  Goldberg  and  Tarjan  [1986]  is  the 
number  of  non-saturating  pushes,  which  dominates  the  number  of  saturating 
pushes  for  the  following  reason:   The  saturating  pushes  cause  structural 
changes  -  they  delete  saturated  arcs  from  the  residual  network.   This 
observation  leads  to  a  bound  of  0(nm)   on  the  number  of  saturating 
pushes  —no  matter  in  which  order  they  are  performed.   The  non-saturating 
pushes  do  not  change  the  structure  of  the  residual  network  and  seem  more 
difficult  to  bound.    Goldberg  [1985]  reduced  the  non-saturating  pushes  to 
0(n-^)  by  examining  the  nodes  in  a  specific  order,  and  Goldberg  and  Tarjan 
[1986]   reduced  the  average  time  per  non-saturating  push  by  using  a 
sophisticated  data  structure.   In  the  next  section  we  show  that  by  using  scaling, 
we  can  dramatically  reduce  the  number  of  non-saturating  pushes  to 
0(n2  log  U) . 
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3.        The  Scaling  Algorithm 

Our  maximum  flow  algorithm  improves  the  generic  preflow-push 
algorithm  of  Section  2,   and  uses  "excess  scaling"  to  reduce  the  number  of 
non-saturating  pushes  from   O(n-m)    to  0{rr  log  U)  .    The  basic  idea  is  to 
push  flow  from  active  nodes  with  sufficiently  large  excesses  to  nodes  with 
sufficiently  small  excesses  while  never  letting  the  excesses  become  too  large. 

The  algorithm  performs  K  =  flog  U1  +  1   scaling  iterations.   For  a 

scaling  iteration,  the  excess-dominator  is  defined  to  be  the  least  integer   A   that 
is  a  power  of  2  and  satisfies     ej  <  A    for  all    i  e  N  .    Further  ,  a  new  scaling 

iteration  is  considered  to  have  begun  whenever    A  decreases  by  a  factor  of  2. 
In  a  scaling  iteration  we  assure  that  each  non-saturating  push  sends  at  least 
A/2   units  of  flow  and  the  excess-dominator  does  not  increase.   To  ensure  that 
each  non-saturating  push  has  value  at  least   A/2  ,  we  consider  only  nodes 
with  excess  more  than  A/2  ;  and  among  these  nodes  with  large  excess,  we 
select  a  node  with  minimum  distance  label  . 

We  show  that  after  at  most   4n-     non-saturating  pushes,  the 
excess-dominator  decreases  by  a  factor  of  at  least  2  and  a  new  scaling  iteration 
begins.   After  at  most  K  scaling  iterations,  all  node  excesses  drop  to  zero  and 
we  obtain  a  maximum  flow. 

In  order  to  select  an  active  node  with  excess  more  than  A/2  and  with  a 

minimum  distance  label  among  such  nodes,    we  maintain  the  lists 

LIST(r)  :  =     (i  e  N  :  ej  >  A/2    and    d(i)  =  r)     for  each  r  =  1,  .  .  .  ,  n-1  .    These 

lists  can  bo  maintained  in  the  form  of  either  a  linked  stack  or  a  linked  queue 
(sec,  for  example,   Aho,  Hopcroft  and  Ullman  [1974]),  which  enables  insertion 


18 


and  deletion  of  elements  in  0(1)  time.)   The  variable  level  indicates  the 
smallest  index  r   for  which   LIST(r)  is   non-empty. 

As  per  Goldberg  and  Tarjan,  we  use  the  following  data  structure  to 
select  efficiently  the  eligible  arc  for  pushing  flow  out  of  a  node.   We  maintain 
with  each  node  i  a  list,  A(i) ,  of  arcs  directed  out  of  it.  Arcs  in  each  list  can  be 
arranged  arbitrarily,  but  the  order  once  decided  remains  unchanged 
throughout  the  algorithm.    Each  node  i  has  a  current  arc  (i,  j)    which  is  the 
current  candidate  for  pushing  flow  out  of  i.   Initially,  the  current  arc  of  node  i 
is  the  first  arc  in  its  arc  list.      This  list  is  examined  sequentially  and  whenever 
the  current  arc  is  found  to  be  ineligible  for  pushing  flow,  the  next  arc  in  the 
arc  list  is  made  the  current  arc.     When  the  arc  list  is  completely  examined, 
we  say  that  the  next  arc  is  null.   At  this  point,  the  node  is  relabeled  and  the 
current  arc  is  again  the  first  arc  in  its  arc  list. 

The  algorithm  can  be  described  formally  as  follows: 
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procedure     MAX-FLOW; 
begin 

PRE-FROCESS; 
K:  =    '  log  U1  ; 
for    k  =  0  to  K  do 
begin 

A  =  2^^-'^ 

for  each    i  e  N     do     if    e,  >  A/2    then  add  i  to  LIST(d(i)); 

level  :  =  1  ; 

while   level  <  2n  do 

if  LIST(level)  =  o    then  level:  =  level  +  1 
else 

begin 

select  a  node  i   from  LIST(level); 
PUSH/RELABEL(i) 
end; 
end; 
end; 
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procedure  PUSH/RELABEL(i); 
begin 

found:  =  false  ; 

let  (i,  j)  be  the  current  arc  of  node  i; 

while  found  =  false    and   (i,  j)  ^  null  do 

if  d(i)  =  d(j)  +  1    and    r-  >  0    then  found:  =  true 

else    replace  the  current  arc  of  node  i  by  the  next  arc  (i,  j); 

if    found  =  true    then 
begin    {push} 

push   min{ej ,  rjj ,  A  -  ej   units  of  flow  on  arc  (i,  j)  ; 

update  residual  capacities  and  the  excesses; 

if  (the  updated  excess)    e^  <  A/2  ,  then  delete  node  i 

from   LIST(d(i)); 
if  j  ;t  s  or  t  and  (the  updated  excess)  e;  >  A/2  ,  then 

add  node  j  to  LIST(d(i))  and  set  level:  =  level  -  1; 
end 


else 


begin  (relabel) 

delete  node  i   from  LIST(d(i)); 

update  d(i):  =  min{d(j)  +  1  ;   (i,  j)  e  A(i)  and  r-  >  0}  ; 

add  node  i  to  LIST(d(i))  and  let  the  current  arc  oi 
the  node  be  the  first  arc  of  A(i); 
end; 


end: 
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4.        Complexity  of  the  Algorithm 

In  this  section  ,  we  show  that  the  distance-directed  preflow-push 
algorithm  with  scaling  correctly  computes  a  maximum  flow  in 
0(nm    +  n-^  log  I"        ne. 

Lemma  5.      The  ithm  satisfies  the  following  two  conditions: 

C3.       Each  non-saturating  push  from  a      node  i    to  a  node  j  sends  at 

least    A/2  units  of  flow. 
C4.       Ko  excess  increases  above    A  .    (i.e.,    the  excess-dominator  does 

not  increase  subsequent  to  a  push.) 

Proof.  For  ever)' push  on  arc  (i,  j)    we  have  ej>A/2    and  e:  <  A/2  ,  since 

node  i   is  a  node  with  smallest  distance  label  among  nodes  whose  excess  is  more 

than  A/2,  and  d(j)  =  d(i)  -  1  <  d(i)  by  the  property  of  the  push  operation.  Hence  , 
by  sending   min  {e^ ,  r-  ,  A  -  e;  )  >  min  {A/2  ,  r-)  units  of  flow,  we  ensure  that 

in  a  non-saturating  push  the  algorithm  sends  at  least   A/2  units  of  flow. 
Further,   the   push  operation  increases  e;    only  .    Let  e':  be  the  excess  at  node  j 
subsequent  to  the  push.    Then   e':  =  e.+    min  {  A/2  ,  r-  ,  A  -  e;  }  <  e.+  A  -  e;  <  A  . 

All  node  excesses  thus  remain  less  than  or  equal  to   A  . 

While  there  are  other  ways  of  ensuring  that  the  algorithm  always 
satisfies   the  properties  stated  in  the  conditions  C3  and  C4,  pushing  flow 
from  a  node  with  excess  greater  than   A/2   and  with  minimum  distance 
among  such  nodes  is  a  simple  and  efficient  way  of  enforcing  these  conditions. 

With  properties   C3  and   C4,   the  push  operation  may  be  viewed  as  a 
kind  of  "restrained  greedy  approach."    Property  C3  ensures  that  the  push 
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from    i   to  j  is  sufficiently  large  to  be  effective.    Propert}'  C4  ensures  that  the 
maximum  infeasibility  of  flow  conservation  never  exceeds    A    during  an 
iteration.   In  particular,  rather  than  greedily  getting  rid  of  all  its  excess,  node  i 
shows  some  restraint  so  as  to  prevent    e:    from  exceeding  A.    (It  is  possible 

that  the  maximum  excess  increases  during  a  push.   In  particular,  if  the 
maximum  excess  is  between  A    and    A/2    ,  it  may  increase  to  A.) 

Keeping  the  maximum  excess  as  low  as  in  C4  may  be  very  useful  in 
practice  as  well  as  in  theory.   Its  major  impact  is  to  "encourage"  flow  excesses 
to  be  distributed  fairly  equally  in  the  network.   This  distribution  of  flows 
should  make  it  easier  for  nodes  to  send  flow  towards  the  sink.    (Alternatively, 
if  a  number  of  nodes  were  to  send  flow  to  a  single  node  j  ,  it  is  likely  that 
node  j   would  not  be  able  to  send  the  accumulated  flow  directly  to  the  sink; 
much  of  the  excess  of  node  j  would  have  to  be  returned  from  where  it  came.) 

Lemma  6.      If  each  push  satisfies  conditions  C3    and  C4,   then  the  number 
of  non-saturating  pushes  per  scaling  iteration  is  at  most    Srr- . 

Proof.  Consider  the  potential  function  F  =   ^  ej  d(i)  . 

i€  N 

The  initial  value  of  F  is  bounded  by  2n-  A    because  ej  is  bounded  by  A  and 
d(i)  is  bounded  by  2n  .    When  the  algorithm  examines  node  i  ,  one  of  the 
following  two  cases  must  apply: 
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Case  I.  The  algorithm  is  unable  to  find  an  arc   along  which  flow  can  be 

pushed.  In  this  case  no  arc  (i,  j)  satisfies      d(i)  =  d(j)  +  1  and  Tj:  >  0  and  the 

distance  label  of  1e  i  goes  up  by    5  >  1  units.    The  increase  in  F  is  thus 

bounded  by    5:  Since  the  initial  value  of  d(i)   and  all  its  increases  are 

bounded  by  2n,  i  le  increases  in  F  due  to  relabelings  of  nodes  are  bounded  by 

2n2A    . 

Case  2.  The  algorithm  is  able  to  identify  an  arc  on  which  flow  can  be 

pushed  and  so  it  performs  either  a  saturating  or  a  non-saturating  push.   In 
either  case,   F  decreases.   Moreover,  each  non-saturating  push   sends  at  least 
A/2  units  of  flow  to  a  node  with  smaller  distance  label  with  a  corresponding 
decrease  in  F  of  at  least  A/2  units.   Since  the  initial  value  of  F  plus  the 
increases  in   F  cannot  exceed  4n-^  A  ,  this  case  cannot  occur  more  than  8n-^ 
times.       (A  more  careful  analysis  leads  to  a  bound  of  4n'^  . )    ■ 

Theorem  1.   The  procedure  Max-Flow  performs   ©(n-^  log  U)   non-saturating 
pushes. 

Proof.      The  initial  value  of  the  excess-dominator  A  is  2'"^°S  ^  "1  <  2U  .    By 
Lemma  5  ,   the  value  of  the  excess-dominator  decreases  by  a  factor  of  2  within 
8n-   non-saturating  pushes  and  a  new  scaling  iteration  begins.   After   1  +  flog 
U  1  such  scaling  iterations,    A  <  1    and  by  the  integrality  of  the  flows  e^  =  0 

for  all    i  e  N  -  {s,  t}.  The  algorithm  thus  obtains  a  feasible  flow,  which  by 
Lemma  4   must  be  a  maximum  flow. 
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Theorem  2.   The  complexity  of  the  maximum  flow  algorithm  is 

0(nm  +  n^  log  U)  . 

Proof:  The  complexity  of  the  algorithm  depends  upon  the  number  of 

executions  of  the  while  loop  in  the  main  program.    In  each  such  execution 
either  a  PUSH/RELABEL(i)  step  is  performed  or  the  value  of  the  variable 
level  increases.     Each  execution  of  the  procedure 
PUSH/RELABEL(i)    results  in  one  of  the  following  outcomes: 

Case  1.  A  push  is  performed.  Since  the  number  of  saturating  pushes  are 
0(nm)  and  the  non-saturating  pushes  are  Oirr  log  U)  ,  this  case 
occurs    0(nm  +  n^  log  U)     times. 

Case  2.    The  distance  label  of  node  i  goes  up.   By  Corollary  2,  this  outcome 
can  occur  0(n)   times  for  each  node  i.     Computing  the  new  distance 
label  requires  examining  arcs  in  A(i).   Since  each  arc  in   A  is 
examined   0(n)   times,  the  total  effort  for  all  of  the  0(n-^)  distance 
increases  is  0(nm)  . 

Thus  the  algorithm  calls  the  procedure   PUSH/RELABEL(i)   0(nm  +  n^  log  U) 
times.   The  effort  needed  to  find  an  arc  to  perform  the  push  operation  is   0(1)   plus 
the  number  of  times  the  current  arc  of  node  i   is  replaced  by  the  next  arc  in  A(i). 
After    I  A(i)  I    such  replacements  for  node  i.  Case  2  occurs  and  distance  label  of  node 

i  goes  up.  Thus,  the  total  effort  needed  is    Zj     2n  lA(i)l    =  0(nm)   times 

iG  N 
plus  the  number  of    PUSIi/RELABEL(i)  operations.   This  is  clearly 
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0(nm  +  n-  log  U).  The  lists  LIST(i)  are  stored  as  linked  stacks  and  queues, 
hence  addition  and  deletion  of  any  element  takes   Oil)  time.     Consequently, 
updating  these  lists  is  not  a  bottleneck  operation.   Finally,   we  need  to  bound 
the  number  of  increases  of  the  variable  level.    In  each  scaling  iteration,   level  is 
bounded  above  by  2n  -  1    and  bounded  below  by  1.  Hence  its  number  of 
increases  per  scaling  iteration  is  bounded  by  the  number  of  decreases  plus  2n  . 
Further,   level  can  decrease  only  when  a  push  is  performed  and  in  such  a  case  it 
decreases  by  1.  Hence  its  increases  over  all  scaling  iterations  are  bounded  by  the 
number  of  pushes  plus  (2n  log  U),   which  is  again   0(nm  +  n-^  log  U)  .    ■ 

5,        Refinements 

As  a  practical  matter,  several  modifications  of  the  algorithm  might 
improve  its  actual  execution  time  without  affecting  its  worst  case  complexity. 
We  characterize  the  modifications  as  one  of  the  following  three  types: 

1.  Modify  the  scale  factor. 

2.  Allow  some  non-saturating  pushes  of  small  amount. 

3.  Try  to  locate  nodes  disconnected  from  the  sink. 

The  algorithm  in  the  present  form   uses  a  scaling  factor  of  2,  i.e.,  it 
reduces  the  excess-dominator  by  a  factor  2   in  the  consecutive  scaling 
iterations.   However,  in  practice  it  might  be  better  to  scale  the  excess- 
dominator  by  some  other  fixed  integer  factor      P  >  2.     The  excess-dominator 
will  be  the  least  power  of  P  that  is  no  less  than  the  excess  flow  at  any  node, 
and  property  C3  becomes 
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C3'.      Each  non-saturating  push  from  a  node  i  to  a  node  j  sends  at  least 
A/p  units  of  flow. 

The  procedure  maxflow  can  easily  be  altered  to  incorporate  the  P  scale 
factor  by  letting  LIST(j)  =  (i:  d(i)    >  A/p  }.  The  algorithm  can  be  shown  to 
run  in  0(nm  +  pn-^  logo  U)    time.   Clearly,  from  the  worst  case  point  of  view 

any  fixed  value  of  P  is  optimum.  The  best  choice  of  the  value  of    P  in 
practice  should  be  determined  empirically. 

Our  algorithm  as  stated  selects  a  node  with  ej  >  A/2    and  performs  a 
saturating  or  a  non-saturating  push.   We  could,  however,  keep  pushing  the 
flow  out  of  this  node  until  either  we  perform  a  non-saturating  push  of  value 
A/2   or  reduce  its  excess  to  zero.   This  variation  might  produce  many 
saturating  pushes  from  the  node  and  even  decrease  its  excess  below  A/2. 

Also,  the  algorithm  as  stated  sends   at  least  A/2  units  of  flow  during 
every  non-saturating  push.    The  same  complexity  of  the  algorithm  is  obtained 
if  for  some  fixed  r  >  1  ,  one  out  of  every  r  >  1   non-saturating  pushes  sends 
at  least  A/2  units  of  flo\v. 

One  potential  bottleneck  in  practice  is  the  number  of  relabels.    In 
particular,  the  algorithm  "recognizes"  that  the  residual  network  contains  no 
path  from  node  i  to  node  t    only  when  d(i)  exceeds  n  -  2  .    Goldberg  [1987] 
suggested  that  it  may  be  desirable  occasionally  to  perform  a  breadth  first 
search  so  as  to  make  the  distance  labels  exact.  He  discovered  that  a  judicious 
use  of  breadth  first  search  could  dramatically  speed  up  the  algorithm. 
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An  alternative  approach  is  to  keep  track  of  the  number   n^    of  nodes 
whose  distance  is  k.     If    n^,    decreases  to  0    after  any  relabel,  then  each  node 
with  distance  greater  than  k   is  disconnected  from  the  sink  in  the  residual 
network.    (Once  node  j   is  disconnected  from  the  sinks,  it  stays  disconnected 
since  the  shortest  path  from  j  to  t    is  nondecreasing  in  length.)   We  would 
avoid  selecting  such  nodes  until  all  nodes  become  disconnected  from  the 
sink.   At  this  time,  the  excesses  of  active  nodes  are  sent  back  to  the  source. 

6.        Conclusions  and  Summary 

Our  improvement  of  the  distance-based  preflow-push  algorithm  has 
several  advantages  over  other  algorithms  for  the  maximum  flow  problem. 

First,  our  algorithm  is  the  most  efficient  algorithm  for  the  maximum 
flo\s'  problem  under  the  reasonable  assumption  that   U  is  polynomially 
bounded  in  n   .   (We  can  relax  this  assumption  on  U    if  we  modify  the 
algorithm  slightly  and  adopt  the  more  realistic  logarithmic  model  of 
computation  as  discussed  in   [19S7a].  ) 

Second,  our  model  has  several  intuitively  appealing  features  that  are 
not  guaranteed  in  other  preflow-push  methods.   We  always  try  to  push  flo^v 
from  a  node  with  a  large  excess.   (Gabow's  scaling  version  of  Dinic's 
algorithm  always  selects  an  augmentation  along  a  path  with  a  large  capacity.) 
This  feature  accounts  for  much  of  the  efficiency  of  his  method.     Also,  we 
never  allow  too  much  excess  to  accumulate  at  a  node.   In  addition,  our 
algorithm  relies  on  verv  simple  data  structures  with  little  computational 
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overhead.    This  feature  contrasts  dramatically  with  the  algorithms  that  rely 
on  dynamic  trees. 

Third,  our  algorithm  is  a  novel  approach  to  combinatorial  scaling 
algorithms.   In  the  previous  scaling  algorithms  developed  by  Edmonds  and 
Karp  [1972]  ,  Rock  [1980]  ,  and  Gabow  [1985]  ,  scaling  involved  a  sequential 
approximation  of  either  the  cost  coefficients  or  the  capacities  and 
right-hand-sides,    (e.g.,  we  would  first  solve  the  problem  with  the  costs 
approximated  by   c/2^  for  some  integer  T.    We  would  then  reoptimize  so  as 
to  solve  the  problem  with  c    approximated  by  c/2'^~^  .    We  would  then 
reoptimize  for  the  problem  with  c  approximated  by  c/2^~^,  and  so  forth.) 
Our  scaling  method  does  not  fit  into  this  standard  framework.    In  Ahuja  and 
Orlin[1987a],  we  describe  an  alternate  framework  for  scaling  algorithms 
which  includes  all  of  the  previous  scaling  algorithms  as  well  as  the  algorithm 
described  here  and  >he  algorithms  of  Goldberg  and  Tarjan  [1987]   for  the 
minimum  cost  flow  problem,  and  the  algorithms  of  Gabow  and  Tarjan  [1987] 
for  the  assignment  and  related  problems. 

Fourth,  our  algorithm  appears  to  be  quite  efficient  in  practice.   In 
Ahuja  and  Orlin  [1987c],   we  empirically  test  a  number  of  different  algorithms 
for  the  maximum  flow  problem.    The  algorithm  presented  here  appears  from 
preliminary  results  to  be  comparable  to  the  best  other  maximum  f\o\v 
algorithm  for  a  range  of  distributions.    In  addition,  our  algorithm  involves 
only  simple  data  structures  so  that  it  is  relatively  easy  to  encode. 
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Fifth,  our  algorithm  can  be  implemented  efficiently  on  a  parallel 
computer  with  a  small  number  of  parallel  processors.   In  particular,  suppose 
that   d  =  m/n  ,  and  that    log  U  =  O(log  n).    Then  our  algorithm  runs  in 
0(n-^  log  n)  time  on  a   parallel  random  access  machine  (PRAM)  using  only 
d  parallel  processors.   This  time  bound  is  comparable  to  the  time  bound 
achieved  by  Goldberg  and  Tarjan  [1986]  using    n  parallel  processors.  Thus  we 
introduce  a   n/d  improvement  in  utilization  of  computational  resources.     In 
addition,  our  model  of  computation  is  less  restrictive.    These  results  are 
discussed  in  Ahuia  and  Orlin  [1987b]. 
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